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Abstract 

We solve the image denoising problem with a dictionary learning technique by writing a convex 
functional of a new form. This functional contains beside the usual sparsity inducing term and 
fidelity term, a new term which induces similarity between overlapping patches in the overlap 
regions. The functional depends on two free regularization parameters: a coefficient multiplying 
the sparsity-inducing L\ norm of the patch basis functions coefficients, and a coefficient multiplying 
the L2 norm of the differences between patches in the overlapping regions. The solution is found by 
applying the iterative proximal gradient descent method with FISTA acceleration. In the case of 
tomography reconstruction we calculate the gradient by applying projection of the solution and its 
error backprojection at each iterative step. We study the quality of the solution, as a function of the 
regularization parameters and noise, on synthetic datas for which the solution is a-priori known. 
We apply the method on experimental data in the case of Differential Phase Tomography. For 
this case we use an original approach which consists in using vectorial patches, each patch having 
two components: one per each gradient component. The resulting algorithm, implemented in the 
ESRF tomography reconstruction code PyHST, results to be robust, efficient, and well adapted to 
strongly reduce the required dose and the number of projections in medical tomography. 



2 



I. INTRODUCTION 



Users of X-ray tomography aim to push the frontiers of their studies towards new domains 
which require finer time resolution, better signal to noise ratio, and less radiation damage. 
All these three requirements bring to a data-starving situation where for a given quality 
goal, the available data volume is never enough. A solution to this problem consists in 
filling in the gap left by the missing data with an a-priori knowledge of the solution. The 
signals occurring in Nature, when cleaned from noise, present most of the time an intrinsic 
sparsity when expressed in the proper basis. An image is intrinsically sparse when it can be 
approximated as a linear combination of a small number n of basis functions with n«iV, 
where N is the image dimensionality. Piece-wise constant images are examples of sparse 
signal : they have non-zero signal only at the flat regions borders when they are expressed 
by their gradient. For piece- wise constant images one can apply very efficient methods based 
on minimization of a convex functional, said also convex objective function, which contains 
a total variation penalty term. For other classes of images such as medical images one has 
to choose different solutions which are adapted to the intrinsic sparsity of the case under 
study. There are mainly two ways : either the sparsity structure is a-priori known and an 
appropriate basis of functions can be built from the beginning, or it must be automatically 
learned from a learning set with the dictionary learning technique. The dictionary learning 
technique has recently been applied to tomography reconstruction using the Orthogonal 
Matching Pursuit (OMP) denoising procedurepQ. This procedure consists in obtaining first 
an over-complete basis of functions and then in fitting every patch of the image using at 
most N omp components selected from this basis. The components are heuristically selected 
choosing, each time, the one having the maximum overlap with the remaining error. This 
optimization method cannot be implemented as a convex objective function optimization 
problem because the linear combination of two candidate solutions can have more than N omp 
components. In other words the optimization domain is not convex. In this paper we present 
a more advanced formalism based on a convex functional that we describe in section [Til For 
the solution of our functional minimization problem we have applied the recently developed 
tools taken from the field of convex optimization|2j. The result is a robust and efficient 



algorithm that we discuss and illustrate with synthetic and medical data in section III 
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II. METHODS 



In this section we introduce first the decomposition of an image into non-overlapping 
patches and the related objective function for denoising. Then we introduce our original for- 
malism which ensures, using overlapping patches, a smooth transition at the patches borders, 
and finally we apply this formalism to image denoising and tomography reconstruction. 

We denote by l p the indicator function of patch p, which is equal to 1 over the patch 
support ( typically an m*m square) and is zero otherwhere. For non-overlapping patches, 
covering the whole domain, we have : 

E ^) = 1 v *- (!) 

p 

where i denotes the pixel position and can be thought as a two-dimensional vector. We 
are looking for the ideal solution x that we express by the vector, w, of its coefficients in 
the basis of patch functions: 

Xi = ^2 5^ w kp ( fk(i ~ r p ). (2) 

p k 

where the set {(fik} is an over-complete basis of functions over the patch support; r p is 
the closest to the origin corner of the patch and Wk p is the component k,p of vector w 
which multiplies the basis function cp^ in the patch p. The denoising problem, given an 
image y, consists in finding the minimum of a functional F(w) = /(w) +#(w) which is sum 
of two terms. The term /(w) = ||y — links the solution to the the data y. The other 
term, g(w), contains the a-priori knowledge about the solution. This way of breaking the 
functional in two terms has his roots in the Bayesian theorem. From a probabilistic point 
of view the denoising problem consists in finding, given a noisy image y of an object, the 
most probable object x that can generate that image. We represent the object x through 
the patches coefficients w. The Bayes theorem, applied to denoising, states that the condi- 
tional probability of w being the exact object given a measurement y, is the product of the 
probability of y being the measure given the exact solution w, times the a-priori probability 
of w. 

Assuming gaussian noise, the conditional probability of y being the measure given the 
exact solution w is exp (— ||y — x^ / (2a 2 )), where x is expressed through the patches co- 
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efficients w by equation [2J The a-priori probability of w is written as exp (—g(w)/ (2a 2 )). 
The a-priori knowledge that the solution is sparse in our patches basis can be expressed by 
using the sparsity-inducing L\ penalization [3]: g(w) = The most probable solution 

w* is obtained by finding the F(w) minimum: 



w* = argmin^(f(w) + #(w)); 

/( w ) = lly - x ll 2 ; #( w ) = z 3 ll w lli • ( 3 ) 

The solution can be obtained by using the iterative shrinkage thresholding algorithm|2j 
(ISTA) iterations: 

w n+ i 7> 7 (w n - 7V/(w n )); w* = w^. (4) 
where T a is the shrinkage operator defined as 

w 

^a(w) Tj — r— max(||w|| 2 — a, 0). (5) 

ll W ll2 

and 7 is a positive number lesser than the inverse of the Lipschitz condition number L: 

7 e]o,i/L]. (6) 

The Lipschitz number L is such that : 

||V/(w 2 ) - V/(wi)|| 2 < L ||w 2 - wi|| 2 ; Vwi, w 2 . (7) 

The ISTA algorithm can be accelerated by the Fast Iterative Shrinkage Thresholding 
gl(FISTA) method. 

In its non-overlapping version, the image denoising with patches is able to detect features 
that are within the field of the patch: if a line crosses the central region of a patch, it will be 
detected if the basis of function has been trained to detect such lines. But in the situation 
where a line intersects only one point in a corner of the patch square, the signal of this point 
is indistinguishable from that of a noisy point, no matters the dictionary training. 

For this reason the patches denoising technique is often used with overlapping patches 
using post-process averaging|5j. In this case the minimization problem is solved for each 
patch separately first, and then the averaging is performed in the overlapped regions. 
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In this study we do not follow this procedure but we add an overlap term into the objective 
function. We choose a system of patches which covers the whole domain, and we allow for 
overlapping. In this case the sum of all indicator functions is greater or equal to one : 

5^1 p (i)>l;Vi. (8) 

P 

We define the core indicator functions 1£, which indicate the core of the patches, and 
make a non-overlapping covering: 

i£(0<ip(0; £ W = i; v.. ( 9 ) 

For a given point z, indicates which patch p has its center C p closest to point z: 

Yj^Wi-cA^Wi-cA^ Vp',i. (10) 

V 

The solution x is composed using the central part of the patches as indicated by the 
functions 1£: 

x i = Yl w kp<Pk(i - r p ). (11) 

p k 

Now we introduce the P operator which is the projection operator, for tomography re- 
construction, and is the identity for image denoising. The functional F(w) whose minimum 
gives the optimal solution is written, for both applications, as: 



F(w) = /(w)+flf(w); ^(w) = /3||w|| 1 ; 
/(w) = ||y-P(x)||? + 

pj^lp(i) I Xi-^2w k p<p k (i-r p ) ) . (12) 

pi \ k / 

where the p factor weights a similarity-inducing term which pushes all the overlapping 
patches, which touch a point z, toward the value Xi of the global solution x in that point. 
For future reference we call ||y — P(x)||2 the fidelity term. 

The solution is found with the FISTA method, using the gradient of /(w) which is easily 
written in compact form: 
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^ = £ 2^(z - r p )l#) <j (P r (P(x) - y)) 4 + 



I Xi - ^2w k , p ,(p k ,(i - ry) ) > + 
v' \ k f J ) 

^2f2(p k {i - r p )pl p (i) I ^2w k > p cp k >(i - r p ) - x { J (13) 

i V k' J 

where P T is the adjoint operator of P, and is called back-projection operator in the case 
of tomography, and is still the identity for image denoising. 



III. APPLICATION 



In this section we compare the fit with patches of a given image with or without overlap- 
ping. We analyze the quality of the recovered image as a function of the two regularization 
parameters f3 and p. The same study is then performed on the same image with additive 
noise. Then we apply the method to tomography reconstruction of virtual phantoms using 
our convex functional with patch overlapping. We analyze the quality of the reconstruc- 
tion as a function of the two regularization parameters projections. Finally, we apply the 
method onto experimental cases for the reduction of number of projection in order to reduce 
the deposited dose during a tomography. 



A. Fitting problem 

Figure [l] shows a training image (Lena) on the left. From this image we have extracted 
by dictionary learning the basis of patches that we show on the right. To obtain the basis 
we have implemented the K-SVD algorithm|6j that we have applied using four atoms in the 
OMP procedure. Our basis consists in an over-complete set of 100 functions having a 7 x 7 
pixels support. 

We apply this set of patches that we have learned from Lena on a completely different 
image (Boy [7j). In figure [2] we show the effect of fitting with non-overlapping patches a 
noiseless image. The original image is subfigure [2^i) . It has been normalized to have its 
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FIG. 1: A training image and its K-SVD basis 

maximum equal to 1 and its minimum equal to 0. From left to right, of the first row, 
in the three columns on the right, we have the results for beta equal to 0.2, 0.02, 0.002. 
The sparsity s is the ratio of non-zero components over the total number of components 
(100 components per every patch here). The value of sparsity for each image is reported 
on its legend. When we take the sparsity average over all the image, it takes the values 
0.029,0.13,0.24 when beta goes from 0.2 to 0.002 in the first row of our example. 

From figure [2] it is obvious how the fit with non overlapping patches gives discontinuity 
at the patches borders. The upper-right corners of the figures are zooms in the hairs zone. 
We can see that, in the case of non-overlapping patches, features in form of lines are more 
sensitive to the regularization term when they are close to the patch border than when they 
are close to the center. When the regularization parameter is weakened by lowering beta, 
then more components are available for the fit and the dimensionality gets high enough to 
fit all the features. 

The second row of figure [2] shows the same fitting problem solved with our overlapping 
patches method. The set of overlapping patches has been obtained by replicating the original 
set of non-overlapping patches by translation vectors (3/,3n); where I and n are integers in 
the interval [0, 2]. This choice of translation vectors results in a nine-fold overlap : each point 
is covered by nine patches. The /3 values take, in the last three columns of the second row, 
the values 0.09,0.009, 0.0009 from left to right. The overlap constraining factor is p = 1. 




FIG. 2: Fitting with no-overlapping (first row) and overlapping (second row) patches a noiseless 
image. The result for f3 = 0.2/9 (e) reproduces the same fidelity/Liratio obtained with/3 = 0.2 (b) 
in the non-overlapping case. For the same level of fitted details the overlapping method requires a 
higher number of components ( bigger s) in order to satisfy the overlapping constraints. 

With these values these images have nearly the same sparsity than the non-overlapping 
results in the same columns of the first row. 

One can see in figure [2] that using overlapping patches the discontinuities have disap- 
peared and the lines are homogeneously represented along their length even for the strongly 
regularized case. The overlap constraint reduces the tessellation effects but we can see that 
for the same sparsity ratio, another kind of distortion, in the form of smoothing, appears. 
This is due to the fact that the fidelity term, which links the data to the solution, is pro- 
portional to the area of the image, while the L\ term acts on all the w p k coefficients whose 
number increases not only with the area of the image but also with the number of replica. 
Taking into account that we have 9 replica per patch, we show in the first image of second 
column the result obtained with f3 = 0.2/9. This rescale back the L\ term to a situation 
which is comparable with the f3 = 0.2 case of image b) of the first row. The level of preserved 
details is comparable in these two cases. The sparsity, however is bigger in the overlapping 
case: s = 0.068 which means nearly seven functions per patch. In the non-overlapping case, 
instead, it was about 3 components per patch. This means that the similarity inducing 
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term is constraining a part of the available degrees to realize smooth transitions between 
neighboring patches. 



B. Denoising problem 

For the denoising problem two things must be kept in mind to select the regularization 
weights: for an image with a stronger noise a stronger value of f3 will be necessary, but it is 
also true that for a stronger value of the regularizing parameter more features of the noiseless 
image will be filtered out. The best value of the regularization parameters is therefore a 
compromise between the necessity of filtering out the noise and that of preserving image 
features. 



a) o=0.1 p=0.016 P=wj^m 






FIG. 3: the effect of overlapping denoising (a,c) on an image with a = 0.1,0.2 (b,d) added white 
gaussian noise. The denoising has been performed with a replication step of 3 and the optimal /3, p 



We define the quality improvement factor Q based on the the Structural SIMilarity (SSIM 
[8]) index S which is 1 when the images are identical. Our definition of Q is : 

Q = (1 - S(y noised) y)) I (1 - S(y denoisedl y)) (14) 

We show in figure [4] the quality improvement factor Q for different values of /?, p for the 
overlapping case and the non-overlapping one. The tests have been done adding a white 
gaussian noise ((J{ no ise} = 0.1 and 0.2) to the original image of figure [2| For the overlapping 
case a replication step of 3 has been used and we have varied separately /3 (squares) and 
p(circles) around the optimal values. In the non-overlapping case (stars) only /3 is varied. 

We observe that the quality improvement factor, for the overlapping case, peaks at a /3 
value which is about 1/9 of the j3 value of the no-overlapping peak. As explained before 
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FIG. 4: The quality improvement factor Q for different values of /3, p around the optimal values for 
° '{noise} =0.1. For the overlapping case a replication step of 3 has been used and we have varied 
separately (3 (squares) and p(circles) around the optimal values. For the non-overlapping case only 
(3 is varied. 

this lower /3 is compensated by the nine-fold increase in the number of components which 
are weighted by Linorm. We note that the optimal quality is better with the overlapping 
method. The trend of /3 as a function of the noise level is as expected : a stronger noise 
needs to be regularized with a stronger f3. 

In figure [3] we show the denoising results (a and c) for two noised images y no i S ed (b 
and d) obtained adding a (7{ no ise} — 0.1,0.2 strong white gaussian to the original image : 
y noised = y original + noise(a{ no i se } ) . The denoising has been performed with a replication 
step of 3 pixels and with the /3 and p giving the optimal SSIM. 
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C. Dose reduction in computed tomography 



In this section we apply the method using overlapping patches for tomography recon- 
struction on synthetic data. The figure [5] presents a 1024 x 1024 pixels phantom and a bases 
of patches that has been learned from the phantom . 

We show in figure [6] the reconstruction obtained from 150 projections between and 180 
degree. Note that this is well below the number of projections, of the order of thousand, 
that would be required in this case by the Shannon-Nyquist sampling theorem. The lack of 
information due to the small number of projections results in a noised reconstruction when 
the filtered-back-projection reconstruction is applied (left) while the image recovered with 
our method (right), with (3/,3n) translations for replication, maintains a high quality. 
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FIG. 5: A phantom and its K-SVD basis 



For this noiseless case the choice for f3 and p was easy and figure [6] was produced with 
our first guess which gave already an excellent reconstruction, without need of further opti- 
mization. The first guess was /3 = 0.001 and p = 50. The quality improvement Q between 
the FBP result and our method is 20. The improvement factor is in this case the FBP error 



(the SSIM index distance from 1 as in 14) divided by the error of our method. 



The noisy case is instead more complex. To study it we add to the phantom sinogram a 
gaussian white noise with a a equal to 5% of the maximum sinogram value. In figure [7] we 
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FIG. 6: reconstruction from 150 projections of a 1024 x 1024 phantom. Top: FBP; Down: our 
method with (3 0.001 and p 50. 

show the Q dependency versus f3 (squares) and p(circles) not far from the optimal choice. 
We have not fully optimized Q because we have found that in this extreme case where p 
tends to be very big, we get a substantial decrease in convergence speed after the shoulder in 
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FIG. 7: the Q dependency versus (3 (squares) and p(circles) not far from the optimal choice of (3 
and p for the reconstruction of the phantom from 150 noisy projection. 



the p (circles) curve. Before the shoulder we have convergence in some hundreds of iterations 
while after the shoulder it takes some thousands to converge. 

We observe that at these high values of p the effective behavior of our method is similar 
to those of a total variation penalization : the image is flattened over large regions. We 
think that the increase in convergence time is due to the fact that the algorithm takes a 
longer time to propagate when the flattened regions are larger. 

Figure [8| bottom, shows the result obtained with our method using the (3 — 1.1 and 
p = 2.5 * 10 5 values obtained from the Q optimization illustrated figure [7J The FBP re- 
construction is reported at the top of the figure. On these images it is clear that that our 
method outperforms standard FBP. The ellipsoids become visible and the plotted profile 
follows the expected behavior, whilst the FBP reconstruction shows only noise. 
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D. Differential Phase tomography 

Finally we show a promising application of our method to medical tomography for a 
sample imaged using X-ray phase contrast imaging(PCI). This technique has shown an 
enhancement of soft tissue visualization in comparison to conventional techniques [9]. PCI 
employs the dual property of X-rays of being simultaneously absorbed and refracted while 
passing through a tissue. Among all the phase contrast techniques we chose to test our 
method on the analyzer based imaging because of the high sensitivity and unique results 
provided by this modality for investigating large and highly absorbing biological tissues (i.e. 
full human breasts) [TU]. Because breast are highly radio sensitive organs, X-ray CT of 
these organs are not clinically applied, even if a 3D would be benefit for radiologists. A 
reduction of the deposited radiation dose in CT combined with the unprecedented contrast 
improvement offered by PCI is thus of high interest for breast cancer detection. 

In the analyzer based PCI technique, the projection data contain a signal which is pro- 
portional to the gradient of the X-ray phase in one direction (i.e. the direction perpendicular 
to the plane formed by the incoming and diffracted X-rays on a perfect Bragg crystal which 
is used for analyzing the radiation passing through the sample). When the object is rotated 
around an axis (Z-axis, for instance), this signal contains contributions from the X and Y 
gradient components, where the X and Y axis corotate with the sample. The two compo- 
nents are de-phased by a rotation angle of 90 degrees and can be reconstructed separately 
by multiplying before-hand the sinogram with the cosine and sine of the rotation angle. We 
apply our formalism considering that the reconstructed and learning images are vectorial 
objects : the value associated with a pixel is not a scalar but a two-component vector. More 
details on the principles and technical aspects of PCI are available in |9j. 

The sample studied is a 7cm human breast imaged with a pixel size of 100 fim. The 
experiment was conducted at the biomedical beam line at the European Synchrotron Radi- 
ation Facility (ESRF). The sample was a human breast mastectomy specimen. The study 
was performed in accordance with the Declaration of Helsinki and was approved by the local 
ethics committee. A monochromatic X-ray beam with energy of 60 keV was used to image 
the breast cancer sample. Result of reconstruction is shown in figure [lO|(top). On this image, 
radiologists could easily identify the skin, fat and glandular tissue. 

The training set is obtained from another breast sample imaged by the same technique 
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but with an high quality reconstruction. We consider a slice image, we apply a Sobel filter 
to extract the two derivative components and we run the KSVD algorithm. 

Figure [9] shows the patches basis functions that we use to fit both components at the 
same time. The patches size is 7 x 7 pixels and each basis function is displayed as a 14x7 
rectangle whose upper 7x7 part is the X component and the lower one the Y component. 



Figure 10 is the reconstruction for a 765 x 765 pixels slice, using only 200 projections 
over the 1000 available. The upper left square is a zoom in the region marked in sub figure 



10, The right column is the reconstruction with our method for X and Y components, 



while the left column is reconstructed with standard filtered back-projection using all 1000 
available projections. Using our method we can still generate a high quality image with 
only one fifth of the projections which would be otherwise necessary to generate a high 
quality reconstruction with the standard FBP method. Visually the difference between 
the FBP results obtained with full data set and our method with a five-fold reduction of 
data is barely noticeable. The different borders of structures like skin layers, fatty tissues, 
and collagen strands are easily identified. The obtained result are very promising and a 
systematic evaluation for clinical application is under-way. The radiation dose absorbed by 
the sample during 200 projections is comparable to that of a standard clinical dual view 
(2D) mammography (3.5mGy). 

For an eventual future clinical application of the PCI method it is important to investigate 
which is the acceptable compromise in terms of low dose and sufficient level of image quality. 
We need therefore to better explore how the quality of the reconstruction is degraded when 
we reduce the dose (i.e. number of projections and the acquisition time) further below the 
standard values. To this end, we performed a reconstruction with only 125 projections and 
results are shown in the figure III-D for one gradient differential image. 

For an eventual future clinical application of the PCI method it is important to investigate 
which is the acceptable compromise in terms of low dose and sufficient level of image quality. 
We need therefore to better explore how the quality of the reconstruction is degraded when 
we reduce the dose (i.e. number of projections and the acquisition time) further below the 
standard values. To this end, we performed a reconstruction with only 125 projections and 
results are shown in the figur^TTj The first column present the result using our method, the 
second column is the result of reconstruction using FBP algorithm. 

If a slightly higher noise level is tolerable, the method may be used with very few pro- 
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jections and thus applied to the screening and diagnosis of human breast cancers with an 
even lower radiation dose than conventional dual mammography. The results of our re- 
construction show an image quality and a capability of discriminating fine structures that 
are still clinically acceptable. On the contrary, images produced with the standard FBP 
reconstruction method are very noisy and not diagnostically satisfactory. 

IV. CONCLUSION 

We have presented a new convex functional which implements in a mathematically pure 
form the concept of overlapping-patches-averaging, which was used so-far with a non-convex 
formalism. The resulting algorithm is robust, efficient, and well adapted to strongly reduce 
the noise in a natural image. The method was applied also to a medical diagnostic case by 
considering phase contrast tomographic data of whole cancer bearing human breasts acquired 
with phase contrast imaging. We demonstrated that it is possible to reduce the deposited 
dose in breast CT by a factor 5 compared to the standard algorithm while keeping the same 
image quality. Although we used this specific example as proof of principle in this study, the 
method we developed and described can be easily applied to other medical tomography fields. 
The numerical results have been generated with PyHST [Hi [12], tthe ESRF tomography 
reconstruction code which uses the GPU implementation of the presented methods. 
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FIG. 8: Reconstruction from 150 angles of the phantom sinogram after adding a gaussian white 
noise with a a equal to 5% of the maximum sinogram value. Top : FBP. Bottom: our method with 
= l.l , p = 2.5* 10 5 . 
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FIG. 9: The vectorial basis of patches learned from a high quality reconstruction of the phase 
gradient of a human breast. In each patch the upper 7x7 part is the X component while the lower 
part is the Y component. 
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FIG. 10: Phase gradient X(first row) and Y(second row) components reconstruction. Left 
with the full set of data. Right: our method with one projection over five, using (5 = 3 * 10 
p=10. 



22 




23 



